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METHOD FOR ANALYTICAL JACOBIAN COMPUTATION IN 

MOLECULAR MODELING 

CROSS-REFERENCES TO RELATED APPLICATIONS 
5 This application is entitled to the benefit of the priority filing dates of 

Provisional Patent Application No. 60/245,730, filed 2000 Nov. 2; and in addition, No. 
60/245,688, filed 2000 Nov. 2; No. 60/245,731, filed 2000 Nov. 2; and No. 60/245,734, filed 
2000 Nov. 2; all of which are hereby incorporated by reference. > 

1 0 BACKGROUND OF THE INVENTION 

The present invention is related to the field of molecular modeling and, more 
particularly, to computer-implemented methods for the dynamic modeling and static analysis 
of large molecules. 

The field of molecular modeling has successfully simulated the motion 
15 (molecular dynamics or (MD)) and determined energy minima or rest states (static analysis) 
of many complex molecular systems by computers. Typical molecular modeling applications 
have included enzyme-ligand docking, molecular diffusion, reaction pathways, phase 
transitions, and protein folding studies. Researchers in the biological sciences and the 
pharmaceutical, polymer, and chemical industries are beginning to use these techniques to 
20 understand the nature of chemical processes in complex molecules and to design new drugs 
and materials accordingly. Naturally, the acceptance of these tools is based on several 
factors, including the accuracy of the results in representing reality, the size and complexity 
of the molecular systems that can be modeled, and the speed by which the solutions are 
obtained. Accuracy of many computations has been compared to experiment and generally 
25 found to be adequate within specified bounds. However, the use of these tools in the prior art 
has required enormous computing power to model molecules or molecular systems of even 
modest size to obtain molecular time histories of sufficient length to be useful. 

There are two sources of computational complexity for molecular modeling 
tasks involving time integration: 
30 1 . The particular molecular model which is used to describe the locations, velocities and 

mass properties of the constituent atoms, the inter-atomic forces between them, and 
the interactions between the atoms and their surrounding environment; and 



-2 - 



# 



2. The particular numerical method used to advance the model through time. Time is 
advanced repeatedly by very short intervals, called timesteps, until a final time has 
been reached. 

Substantial work has been completed in reducing the computational load for 
5 molecular models, such as the reduction of model complexity by constraining higher order 
modes with rigid body assumptions, simplifying the model with rigid or flexible 
substructuring, Order(/V) dynamics, efficient implicit solvent models, and multipole methods 
for the force field models (see, for example, U.S. Patent No. 5,424,963 on the commercial 
MBO(N)D software package). Co-pending applications, U.S. Appln. No. , 

10 entitled "METHOD FOR LARGE TIMESTEPS IN MOLECULAR MODELING," and U.S. 

Appln. No. , entitled "METHOD FOR RESIDUAL FORM IN MOLECULAR 

MODELING," both of which are filed of even date, claim priority from the previously cited 
provisional patent applications and which are assigned to the present assignee, and which are 
incorporated by reference herein, describe further improvements in molecular models and 

15 numerical methods. 

The primary reason molecular simulations are so slow is that current 
numerical methods require very small timesteps, typically between 1 and 10 femtoseconds 
(10" 15 to 10" 14 seconds). Each timestep taken requires the computation of a new state 
(position and motion for each atom) of the particular molecular model, and then computation 

20 of the new set of forces resulting from the new state. For example, molecular dynamics 
simulations of the complex behavior of large molecules, such as the folding of a protein, 
typically need to cover a time span from at least a microsecond up to several seconds or even 
minutes. With techniques currently in common use, this results in the requirement to take 10 9 
to 10 16 timesteps in the computer simulation. The per-step computations of the state, and 

25 especially the forces, grow very expensive as the problem size increases. Even with the 
fastest computers available today, months, years or even centuries of computer time are 
required to solve such problems even for systems of modest size. 

Heretofore, it has been widely believed by molecular dynamicists that these 
small timesteps are an inevitable requirement of the need to maintain accuracy in the 

30 presence of the very high frequencies to be found in vibrations of molecular bonds. For 
example, see Leach, Molecular Modelling Principles and Applications, 1996, p. 328; * 
Berendsen, in Computational Molecular Dynamics: Challenges, Methods, Ideas Deuflhard et 
al. (ed.), Springer, 1999, pg. 18; Rapaport, The Art of Molecular Dynamics Simulation, 
Cambridge, 1995, reprinted with corrections 1998, p. 57; and U.S. Patent No. 5,424,963. 
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This common-sense belief is incorrect, however. The computer science sub- 



discipline of numerical analysis has produced an extensive theory of numerical integration for 
problems in which high frequencies exist but are of little interest. These problems are termed 
"stiff problems (see, for example, Hairer and Wanner, Solving Ordinary Differential 
5 Equations II: Stiff and Differential-Algebraic Problems, 2 nd ed., Springer, 1996). In these 
cases, it is the stability of the integration method, not the required solution accuracy, which 
limits the timestep. Integration methods exist which have unconditional stability, meaning 
that in theory they can take arbitrarily large timesteps. Such methods have a mathematical 
property called "L-stability." Only so-called "implicit" integration methods can be L-stable, 
10 but very few implicit integration methods actually are L-stable. Previously cited co-pending 



MOLECULAR MODELING," covers the use of implicit and in particular L-stable 
integration methods. 

The present invention covers another critical aspect of allowing the implicit 

15 and in particular L-stable integrators to take large timesteps, namely, more accurate methods 
for computing required components of the implicit integration methods called "Jacobians". 

But the same problem of high frequency molecular vibration for MD 
simulations causes problems for the calculation of Jacobians. For example, the repulsive van 
der Waal's forces that are generated as the electron orbitals of two atoms begin to overlap 

20 must be accounted for in a molecule. This quantum mechanical effect is often approximated 
by the so-called Lennard- Jones potential (Rapaport, op. cit.), which models the repulsive 
force as being proportional to 1/r 13 , where r is the distance between the centers of the 
atoms. This is an extreme stiff interatomic force which is characteristic of molecular 
dynamics (MD) simulations and poses particular difficulty for any numerical integration 

25 scheme used to advance time during the simulation. As a result and as mentioned previously, 
most prior art MD integration schemes have timesteps limited by the high frequencies 
generated by these extremely stiff repulsive forces. And in particular, the stiffness of the 
atomic forces greatly magnify the difficulty of forming certain required Jacobian matrices. 
Such Jacobian matrices are a necessary ingredient of any stable implicit integration scheme, 

30 such as described in the immediately cited co-pending application. 

All general-purpose simulation codes provide routines to compute Jacobians 
numerically as follows. For a given equation to integrate y = f{y,t), the desired Jacobian is 
J = df I dy and is numerically computed: 



U.S. Patent Appln. No. 



, entitled "METHOD FOR LARGE TIMESTEPS IN 
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Ay 

where 

AT = / U -/ U, 

Note that the perturbation Ay has to be selected to give a good result and may 
5 be difficult to choose, especially for stiff systems. In contrast, analytic Jacobians are 

computed by solving directly, or in this case algorithmically, for the equation of the desired 
derivatives. 

Linearized models are regularly produced analytically for simple systems. 
Such linearization is usually performed around an equilibrium solution. It is common in such 

10 packages as ACSL (Advanced Continuous Simulation Language), (ACSL Reference Guide, 
Mitchell Gauthier and Associates, 1996), and SPICE (a circuit simulation package), (R. 
Kielkowski, Inside SPICE, McGraw-Hill, 1998) to perform equilibrium analysis followed by 
linearization. Such linearization is performed numerically. 

In contrast, the Jacobian of the present invention represents linearization about 

15 an instantaneous solution of the differential equations (non-equilibrium) and is generated 
analytically. It should be noted that another prior art approach to generating analytic 
Jacobians is to use automated differentiation tools such as ADIFOR (Automatic 
Differentiation of Fortran) (C. Bischof, et. al., ADIFOR 2.0 Users ' Guide, Argonne National 
Laboratory, 1998) that can symbolically differentiate arbitrary equations. These tools could 

20 be used to implement this invention in practice. However, the structure of the equations must 
be exploited properly to minimize the computations, to avoid errors due to round off and term 
cancellations, and to avoid "equation swell" which could limit the size of problems solved. 

Rather, the present invention allows for the calculation of analytic Jacobians 
for the effective implicit integration, including L-stable integrators, of the equations of 

25 motion of molecular models. 

SUMMARY OF THE INVENTION 
The present invention provides for a method of modeling the behavior of a 
molecule. The method has the steps of selecting a torsion angle, rigid multibody model for 
30 said molecule, the model having equations of motion; selecting an implicit integrator; and 

generating an analytic Jacobian for the implicit integrator to integrate the equations of motion 
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so as to obtain calculations of the behavior of the molecule. The analytic Jacobian J is 
derived from an analytic Jacobian of the Residual Form of the equations of motion and is 
described as: 




; and 



J, 



dq _d(Wu) 




dq dq 




where q are the generalized coordinates, u are the generalized speeds, W is a 



joint map matrix and Mis the mass matrix and p u is the dynamic residual of the equations of 
motion, and z is -M~ ] p u (q, u, 0) . The method can also be used for any physical system 

which can be modeled by a torsion angle, rigid multibody system. 

The present invention also provides for the computer code for simulating the 
behavior of a molecule, or any physical system, which can be modeled by a torsion angle, 
rigid multibody system. A module in the computer code with an implicit integrator includes 
the analytic Jacobian as described above. 



Fig. 1 is a representational block module diagram of the software system 
architecture in accordance with the present invention; 

Fig. 2 illustrates the tree structure of the multibody system of the molecular 
model according to the present invention; 

Fig. 3 illustrates the reference configuration of the Fig. 2 multibody system; 

Fig. 4A illustrate a sliding joint between two bodies of the Fig. 2 multibody 
system; Fig. 4B illustrate a pin joint between two bodies of the Fig. 2 multibody system; Fig. 
4C illustrate a ball joint between two bodies of the Fig. 2 multibody system; 

Fig. 5 summarizes general computational steps for the Residual Form method 
and Direct Form methods used for the analytic Jacobian computation; 



BRIEF DESCRIPTION OF THE DRAWINGS 



Fig. 6 is a chart which summarizes the general computational steps for the 
analytic Jacobian; 

Fig. 7 is a plot of digits of accuracy versus perturbation to show the accuracy 
of analytic Jacobian over the numerical Jacobian. 

DESCRIPTION OF THE SPECIFIC EMBODIMENTS 

GENERAL DESCRIPTION OF THE PRESENT INVENTION 

The numerical method used to advance time in the simulation of a modeled 
molecular system is called the integrator, or integration method. The integration proceeds by 
discretizing the governing equations of motion of the molecular system. In the case of an 
implicit integrator, this step results in a set of coupled nonlinear algebraic equations (the 
"stage" equations) and the solution of these equations yields the state vector for the molecular 
system at the next time step. 

The present invention aids the solution of the stage equations. Because the 
atomic forces vary so rapidly over short distances, it is difficult to propagate atomic 
trajectories accurately. Small errors in the atomic positions lead to gross errors in the 
satisfaction of the stage equations. Since the Jacobian is used solve the stage equations 
iteratively, an inaccurate Jacobian leads to trial solution that are far from the correct solution. 
This leads to retraction of trial solutions and hinders the simulation. 

Numerical Jacobians may be correct in only half their significant digits. An " 
analytical Jacobian will often be correct in all but the last digit. The benefit of this result is 
that the integrator can take far fewer time steps to simulate the specified interval, allowing 
full exploitation of the theoretical stability of the integration method. 

The ease or difficulty in producing the Jacobian depends crucially upon the 
formulation used to produce the governing equations. For instance, global Cartesian 
formulations produce equations with very limited explicit coordinate dependence. Producing 
an analytic Jacobian for such a formulation is well known. On the other hand, using internal 
coordinates (in which each molecular subunit's location is expressed relative to an earlier 
subunit's location) as independent variables has great benefits. This method is most valuable 
for any molecule modeled with any choice of internal coordinates, and in particular when 
used with protein models or other polymeric molecules using "torsion angles" between the 
residues as internal coordinates. An algorithm for producing an analytic Jacobian for a 



system formulated in torsion angles is extremely difficult to devise. However, the present 
invention achieves this task. 

The present invention addresses a seemingly intractable problem: production 
of the analytic Jacobian for a formulation using internal coordinates, and specifically torsion 
5 angles, which is generally thought to be impractical. In addition to being more accurate than 
numerical Jacobians, the analytic Jacobians are also cheaper (in computing power) to 
produce. The present invention also recognizes a key result that the Jacobian of the state 
derivatives can be computed by applying a matrix inverse to the Jacobian of the computed 
torque method. This result allows significant savings in computer time and effort to construct 
10 the algorithm. 

Furthermore, this method of producing analytic Jacobians for multibody 
system formulations using torsion angle, internal coordinates has not been seen in the general 
MBS literature. The present invention can be used for any torsion angle MBS formulation, 
which can be applied to many other disciplines besides molecular simulations, including, but 
15 not necessarily limited to, mechanical systems, robotic systems, vehicle systems, or any other 
system which could be described as a set of hinge-connected rigid bodies. 

OVERVIEW OF DESCRIPTION 

The preferred embodiment is divided into several sections. The first set of 

sections describes the MD simulation architecture, the multibody system (MBS) definitions, 
20 and Residual Form of the MBS equations for the subsequent descriptions. The next set of 

sections discusses the definition of the Jacobian, its role in the implicit integration method, 

and the computation of the analytic Jacobian using the Residual Form. Also shown is the 

superior accuracy and performance of the analytic Jacobian vs. the numerical Jacobian. 

Further efficiencies in the computation of the analytic Jacobian are discussed, specifically, 
25 exploiting the rigid body MBS to "contract" the size of the Jacobian matrix, and exploiting 

the topological structure of the MBS to eliminate unnecessary computations. 

To solve ordinary differential equations (ODEs), most of the prior art have 

used equations expressed in the Direct Form, i.e., y = f(yj) (where v means — ). The 

dt 

equations of motion for a biomolecular system can be cast into this form (and called the 
30 Direct Form). In molecular modeling, all prior art known to the present inventors have used 
the Direct Form. That is, q = Wu , u = M' ] f , where q and u are generalized coordinates and 
speeds respectively, so that conventional ODE solution methods can be applied. However, 
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this requires a matrix inversion of M (representing the mass of the system) at a cost of 
Order( TV ) to Order( N 3 ) floating point operations (depending on algorithm used, where N is 
the number of degrees of freedom in the system), since the natural form of the equations 
gives rise to inertial coupling between the derivatives of the generalized speeds. That is, the 
5 equations are most naturally produced in the form q - Wu = 0 , Mil - / = 0 , where M , the 
mass matrix, depends explicitly upon the generalized coordinates i.e., M = M (q) . This 

fact requires forming and effectively factoring the mass matrix each time the state derivatives 
are needed by the integrator in integrating the equations of motion over time. The 
generalized joint map matrix W is block diagonal and, although it is also dependent on the 
10 coordinates W = W (#) , it does not have a significant computational cost. 

In accordance with the present invention, a method for the solution of the 
equations of a molecular system is expressed in Residual Form to bypass the customary step 
of producing the state derivatives directly. The Residual Form method has the following 
steps: 

15 1) Discretization of the solution variables. The specific form of discretization is dictated 

by the particular implicit integration method used to advance the molecular model in 
time. Implicit integration follows from the Residual Form. Implicit integration, 
especially L-stable integrators and other highly stable integrators, such as implicit 
Euler, Radau5, SDIRK3, SDIRK4, other implicit Runge-Kutta methods, and DASSL 

20 or other implicit multistep methods, also provide other advantages for molecular 

modeling. See, for example, the above-cited U.S. Patent Appln. No. , entitled 

"METHOD FOR LARGE TIMESTEPS IN MOLECULAR MODELING," filed of 
even date. As a particularly simple example, when used with implicit Euler 
integration, the discretization is as follows: 

25 q = { c in-q„-i)ih, ii = {u n -u n _ x )/h 

where h is the timestep. 
2) Substitution into the residual equations: 

f P q \ ( q-W{q)u ^ 



M(q)u-f(t 9 q 9 u)J 



3) Solution of the resulting nonlinear algebraic equations 



(p \ 



= 0 for q n and u n 
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The kinematic residual p q compares an estimated q generated from the 

implicit integrator to the derivatives computed by the routines for determining the joints of 
the molecular model, which is described in greater detail below. The second row of the 
residual is p u , the dynamic residual, which determines the degree to which an estimated it 

5 satisfies the equations of motion. 

The system mass matrix M and the so-called 'bias-free hinge torque' / are 
both state dependent. The bias-free hinge torque is generated by the dynamic residual routine 
when the calculated it vector passed to the residual routine is zero. In general, the hinge 
accelerations are a response to applied forces, joint torques, and motion-induced effects (such 
10 as Coriolis and centrifugal forces.) If the system were at rest, and subjected only to joint 

torques, it would be considered in a bias-free state. The real system with its actual inputs can 
be reduced to a bias-free state by computing a set of joint torques equivalent to the biased 
inputs. Both sets produce the same hinge accelerations. 

The preferred embodiment of the Residual Form is shown for an Order( N ) 
15 torsion-angle, rigid-body for the molecular model. The following sections develop the 
molecular model from basic definitions and show how the model is used to compute the 
motion of the model. First, the overall computer code architecture for the molecular model 
simulation is described. Then an Order( TV ) torsion- angle, rigid multibody system is derived, 
along with notation used, the reference configuration, the definitions of the joints between the 
20 bodies, generalized coordinates, and generalized speeds. This approach for dynamics is 
similar to that used by T. R. Kane (Dynamics, 3 rd ed., 1978.) 

Finally, the preferred embodiment of the analytic Jacobians using the results 
of the Residual Form is developed. 

MOLECULAR DYNAMICS SIMULATION ARCHITECTURE 
25 The general system architecture 48 of the software and some of its processes 

for modeling molecules in accordance with the present invention are illustrated in Fig. 1 . 

Each large rectangular block represents a software module and arrows represents information 

which passes between the software modules. The software system architecture has a modeler 

module 50, a biochem components module 52, a physical model module 54, an analysis 
30 module 56 and a visualization module 58. The details of some of these modules are 

described below; other modules are available to the public. 

The modeler module 50 provides an interface for the user to enter the physical 

parameters which define a particular molecular system. The interface may have a graphical 
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or data file input (or both). The biochem components module 52 translates the modeler input 
for a particular mathematical model of the molecular system and is divided into translation 
submodules 60, 62 and 64 for mathematical modeling the molecule(s), the force fields and 
the solvent respectively of the system being modeled. There are several modeler and 
5 biochem components modules available including, for example, Tinker (Jay Ponder, 

TINKER User's Guide , Version 3.8, October 2000, Washington University, St. Louis, MO). 

With the translated physical parameters from the biochem components module 
52, the physical model module 54 defines the molecular system mathematically. At the core 
of the module 54 is a multibody system submodule 66. The analysis module 56, which 

10 communicates with the physical model module 54 and the visualization module 58, provides 
solutions to the computational models of the molecular systems defined by the physical 
model module 54. The analysis module 56 consists of a set of integrator submodules 68 
which integrate the differential equations of the physical model module 54. The integrator 
submodules 68 advance the molecular system through time and also provide for static 

15 analyses used in determining the minimum energy configuration of the molecular system. 
The analysis module 56 and its integrator submodules 68 contain most of the subject matter 
of the present invention and are described in detail below. 

The visualization module 58 receives input information from the biochem 
components module 52 and the analysis module 56 to provide the user with a three- 

20 dimensional graphical representation of the molecular system and the solutions obtained for 
the molecular system. Many visualization modules are presently available, an example being 
VMD (A. Dalke, et aL, VMD User's Guide . Version 1.5, June 2000, Theoretical Biophysics 
Group, University of Illinois, Urbana, Illinois). 

The described software code is run on conventional personal computers, such 

25 as PCs with Pentium III or Pentium IV microprocessors manufactured by Intel Corporation of 
Santa Clara, California. This contrasts with many current efforts in molecular modeling 
which use supercomputers to perform calculations. Of course, further speed improvements 
can be obtained by running the described software on faster computers. 

MOLECULAR MODEL AND MULTIBODY SYSTEM DESCRIPTION 
30 The integrators described below in the submodule 68 operate upon a set of 

equations which describe the motion of the molecular model in terms of a multibody system 
(MBS). To aid the computation of the integration methods described in detail below, a 
torsion angle, rigid body model is used to describe the subject molecule system, in 
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accordance with the present invention. Internal coordinates (selected generalized coordinates 
and speeds) are used to describe the states of the molecule. 

The MBS is an abstraction of the atoms and effectively rigid bonds that make 
up the molecular system being modeled and is selected to simplify the actual physical system, 
5 the molecule in its environment, without losing the features important to the problem being 
addressed by the simulation. With respect to the general system architecture illustrated in 
Fig. 1, the MBS does not include the electrostatic charge or other energetic interactions 
between atoms nor the model of the solvent in which the molecules are immersed. The force 
fields are modeled in the submodule 62 and the solvent in the submodule 64 in the biochem 

10 components module 52. 

Fig. 2 illustrates the tree structure of the MBS of a subject molecule. The 
basic abstraction of the MBS is that of one or more collections of hinge-connected rigid 
bodies 1 70. A rigid body is a mathematical abstraction of a physical body in which all the 
particles making up the body have fixed positions relative to each other. No flexing or other 

1 5 relative motion is allowed. A hinge connection is a mathematical abstraction that defines the 
allowable relative motion between two rigid bodies. Examples of these rigid bodies and hinge 
connections are described below. 

One or more of the bodies, called base bodies 172, have special status in that 
their kinematics are referenced directly to a reference point on ground 174. The system graph 

20 is one or more "trees". An important property of a tree is that the path from any body to any 
other body is unique, i.e., the graph contains no loops. The bodies in the tree are n in number 
(the base has the label 1). The bodies in the tree are assigned a regular labeling, which means 
that the body labels never decrease on any path from the base body to any leaf body 176. A 
leaf body is one that is connected to only a single other body. A regular labeling can be 

25 achieved by assigning the label n to one of the leaf bodies 178 (there must be at least one). If 
this body is removed from the graph, the tree now has n - 1 bodies. The label n - 1 is then 
assigned to one of its leaf bodies 180, and the process is repeated until all the bodies have 
been labeled. This is also done for any remaining trees in the system. 

To help maintain the relationship between the bodies, an integer function is 

30 used to record the inboard body for each body of the system. The inboard body for each base 
is ground and /, the parent or inboard body 182 for body k 184, is referred to as i = inb(k) . 
Additionally, the symbol N refers to the inertial, or ground frame 174. A superscript O refers 
to the ground origin (0,0,0). 
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The symbol for the vector from one point to another contains the name of the 
two points. Thus, r PQ is the vector from the point P to point Q. A vector representing the 
velocity of a point in a reference frame contains the name of the point and the reference 
frame: "v^ . Certain symbols to be introduced later relate two reference frames. In this case, 
5 the symbol contains the name of two frames. Thus, i C k is the direction cosine matrix for the 
orientation of frame k in frame i. This symbol refers to the direction cosine matrix for a 
typical body in its parent frame. Thus, l C k (j) indicates the actual body j in question. The 
left and right superscripts do not change with the body index. This is also true for the other 
symbols. An asterisk indicates the transpose: H\k) , for example. A tilde over a vector 
10 indicates a 3 by 3 skew-symmetric cross product matrix: vw = vxw. is an i by i identity 
matrix., and 0,. is a zero vector of length i and £. is an i by i zero matrix. 

Rigid Bodies of the Model 

Fig. 3 illustrates the reference configuration 190 of a sample "tree" of the 
MBS. More than one tree is allowed. A point of each body is designated as Q, its hinge 

15 point. For example point Q k 186 is the hinge point for body k 184. A fixed set of coordinate 
axes is established in the inertial frame 198. An arbitrary configuration of the MBS is chosen 
as its reference configuration 1 90. While in this configuration the image of the inertial 
coordinate axes is used to establish a set of body-fixed axes in each body. In the reference 
configuration each hinge point Q is coincident with P, a point of its parent body (or extended 

20 body.) For each body, point P is called the body's inboard hinge point. So the inboard hinge 
point P k 188 for body k 184 is a point fixed in its parent body i 182. The inboard hinge point 
for each base body is a point O 192 fixed in ground. The expanded view that shown in Fig. 2 
more clearly shows that point Q k 186 is fixed in body k 184 and point P k 188 is fixed in 
parent body i 182. 

25 The hinge point locations define d(k) 194, a constant vector for each body, 

and can also be written r Qi?k . The vector for body k is fixed in its parent body i. It spans 
from the hinge point for body i to the inboard hinge point for body k. The vector d(l) 196 
spans from the inertial origin to the first base body's inboard hinge point (also a point fixed in 
ground), and can be written r 0Q " . 
30 For a body, m(k) , p(k) , and l^ k {k) define the mass properties of body k for 

its hinge point Q k . These are, respectively, the mass, the first mass moment, and the inertia 
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M(k) = 



matrix of the body for its hinge point in the coordinate frame of the body. For a rigid body 
made up of a distribution of particles, the mass properties are constants that are computed by 
a preprocessing module. The details of these computations can be found in standard 
references, such as Kane, T.R., Dynamics, 3 rd Ed., January 1978, Stanford University, 
Stanford, CA. 

Let M(k) , the spatial inertia of body k for its hinge point Q k , be given by the 
symmetric 6 by 6 matrix 

P(*) ' 
-p(/c) rn(k)£ 3 

Each joint in the system is described by geometric data. For instance, a pin 
10 joint is characterized by an axis fixed in the two bodies connected by the joint. The particular 
data for a joint depends on its type. The number n, the inb function, the system mass 
properties, the vectors d(£), and the joint geometric data (including joint type) constitute the 
system parameters. 

Joints and Generalized Coordinates of the Model 

15 Figs. 4A-4C illustrate the joint definitions of the preferred embodiment of the 

MBS: the slider joint 100, the pin joint 102, and the ball joint 104. Each joint allows 
translational or rotational displacement of the hinge point Q k 106 relative to the inboard hinge 
point P k 108. These displacements are parameterized by q(k) 1 10, the generalized 
coordinates for body k. In passing, it should be noted that generalized coordinates are 

20 examples of generalized quantities, which refer to quantities that have both rotational 
character and translational character. For instance, a generalized force acting at a point 
consists of both a force vector and a torque vector. The generalized coordinate q(k) for the 
slider joint 100 is the sliding displacement x 112. The generalized coordinate q(k) for the 
pin joint 102 is the angular displacement O 114. The generalized coordinate q(k) for the 

25 ball joint 104 is the Euler parameters {s }9 s 29 s^ 9 s 4 ) 116. 

Each joint may be a pin, slider, or ball joint; or a combination of these joints. 
Many other joint types are possible, including, but not limited to, free joints, U-joints, 
cylindrical joints, and bearing joints. For instance, q(k) = (x,^,z), the inertial measure 

numbers of the vector from the base body inboard hinge point to the base body hinge point 
30 express the base body displacement in ground as three orthogonal slider joints. A free joint 
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consists of three orthogonal slider joints combined with a ball joint, and has the full 6 degrees 
of freedom. 

The collection of generalized coordinates for all the bodies comprises the 
vector q , the generalized coordinates for the system. 

5 Given the generalized coordinates for a particular joint, two quantities: 

r PkQk (k) , the joint translation vector and l C k {k) , the direction cosine matrix for body k in its 
parent are formed. The translation vector r PkQk (k) expresses the vector from the inboard 
hinge point P of body k to the hinge point Q of body k, in the coordinate frame of the parent 
body. Details of these computations depend on the joint type and can be easily derived. For 

10 purposes of this description, access to a function that can generate r PkQk (k) and l C k (k) given 
the system generalized coordinates is assumed. 

As introduced, the choice of hinge point for each body is arbitrary. However, 
judicious choice greatly simplifies matters. For instance, for pin joints the hinge point should 
be chosen as a point on the axis of the joint. For this choice points P and Q remain coincident 

15 for all values of the joint angle, so the joint translation is zero. If the point Q is chosen at a 
distance from the axis, points P and Q move relative to each other: 

r PkQk (k) = kx r OQk sin 0 - (1 - cos 0) (f 3 - iX ) r OQk 
where k is the joint axis unit vector, 0 is the joint angle, and r OQk is the vector from any point 
on the axis to point Q. 

20 For pin joints and ball joints, a point on the axis is always chosen as the hinge 

point. For these joints the translation vector r PkQk (k) is zero. 

For a slider joint, the translation vector r PkQk (k) is q(k)*k . 
The direction cosine matrix for a pin is 

'C* (k) = E 3 cos0+k sin 0 + JJt* (1 - cos 0) . 
25 The direction cosine matrix for a slider is £ 3 . 

Generalized Speeds of the Model 

Let l V k (k) , the generalized velocity of the hinge point of body k measured in 
its parent /, be parameterized by u(k) , a set of generalized speeds. Then: 
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H(k) = 



Here, the matrix H(k) is called the joint map for this joint. It is a n u (k) by 6 matrix, where 
n u (k) is the number of degrees of freedom for the joint (1 for a pin or slider, 3 for a ball, 6 for 
a free joint). H(k) can, in general have dependence on coordinates q . Given the 
generalized speeds for the joint, the joint map generates the joint linear and angular velocity, 
expressed in the child body frame. The following are used for the joints: 

H(k) = [A 0 0 0], pin 
H(k) = [0 0 0 X\ 9 slider 
H(k) = [g 3 £3], ball 

£3 i C k (k)_ 

The collection of generalized speeds for all the bodies comprises the vector w, 
the generalized coordinates for the system. As before, access to a function that can generate 
the vector l V k (k) given (q,u) and a specific joint type, is assumed. Access to a function that 
can compute the derivatives q(k) = q(q(k),u(k)) is also assumed. This routine generates the 
time derivative of the generalized position coordinates: 

q = W(q)u 

where W{q) is a block diagonal matrix that relates q and u , with each block depending 
upon the joint type: 

q = u for pin joint, slider joint 



, free 



£4 ^3 <^2 

— £ 2 £ x S 4 
-£ x -£ 2 £ 



CO, 



CO, 



CO, 



for ball joint 



where q = [s } £ 2 £ 3 £ 4 ] and u = [co { co 2 co 3 ]* 



20 



and a free joint is a combination of 3 slider joints and one ball joint. Note that there are 
(derivatives of the Euler parameters) associated with 3 u 's for ball joints. 

Similarly, 1 A k (k) , the generalized acceleration of the hinge point of body k in 
its parent, is given by: 



'A\k) = 



= H'(k)u(k) 



It is these generalized coordinates q, and generalized speeds «, the internal coordinates for 
purposes of this description, of the molecular system which are calculated. Rather than 
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working with the typical inertial coordinates (x , y, z) and speeds in these inertial coordinate 
systems, calculations for the subject molecular system are reduced. 

CALCULATIONS OF THE EQUATIONS OF MOTION 

With the exemplary rigid multibody, torsion angle model described, the 
equations of motion can now be calculated. In accordance with the present invention, the 
motion of the MBS molecular model is determined by the Residual Form. The Residual 
Form method requires calculations termed the "first" kinematic calculations to distinguish 
them from the "second" kinematic calculations, which are further required by the Direct 
Form (which is included in this description for purposes of comparison). 

First Kinematic Calculations for the Molecular Model 

In the first kinematic calculations, given the internal coordinates of the 
molecular system, (<?, w,w) and the system parameters, the following position, velocity and 
acceleration kinematics are computed for each rigid body k of the molecular model. (In 
passing, it should be noted that when the First Kinematic calculations are done for the 
Residual Form method, the it is passed in as a guess of the solution which the integration 
method then refines to the correct solution. In contrast, u is set to zero when used for the 
Direct Form method. This is shown clearly in the later descriptions of the two methods.) 

For each body k compute: 

"C k (k),r Q &(k),r OQ >(k)j0 k (k), 
N o k {k\ N v Qk {k\V{k\ 
N a k (k), N a Q > (k) 9 A(k) 

These computations are done recursively, starting from each base body and progressing to the 
leaves. 

N C k {k) , the direction cosine matrix for body k in ground is defined as: 
iV C*(l) = 'C*(l) 

»C k (k) = N C k (iyC k (k), k = 2 9 ...n, i = inb(k) 
l C k (k) comes from the joint routine described above. 

r QiQk (k) , the position vector from Q. , the hinge point of the parent of body k 

to Q k , the hinge point of body /c, expressed in the parent frame, is defined as: 

r PkQk (k) comes from the joint routine. 
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r OQk (k) , the position vector from the inertial origin O to Q k , the hinge point 
of body k, expressed in the global frame, is defined 
A-°&(l) = r aa (l) 

r° a (k) = r OQk (i) + N C" (/> aa (k), k = 2,...n, i = inb(k) 
'0 k (k) , the rigid body transformation operator for body k is defined 



n c k (k) p Q& (kyc k (ky 

'C k (k) j 

V{k) , the spatial velocity for body k at its hinge point, expressed in the frame 
of body k, is defined 



'#*(*) = 



k = 



V{1) = 



% a (l)J 



= 'V k {\) 



^(A:) , the spatial acceleration for body k at its hinge point, expressed in the 
frame of body k, is defined 



^4(1) = | 
A(k)± 



r N a k {\) 



= '^(1) 



V 



= /4 + 



2a> 



'>*(£) + /t = 2,...«, i = inb(k) 



where 



( 



A = i 0 k '(k)A(i) 



O3 



i C k \k)( N co k (i) x x r aa 

0 = 'C**(*) A '^(i) 

Of course, the computations can all be computed in a single pass if desired. 

After completing these steps for one incremental time step, the MBS can 
service kinematics requests to compute the (generalized) position, velocity, or acceleration 
information for any point of any body. This is done by computing the required information 
for any point in terms of the hinge quantities for its body, using standard rigid body formulas. 

Residual Computation 

With the first kinematic calculations described above, the residual 
computation for the Residual Form method can be determined. A detailed description of the 
Residual Form and its application to molecular modeling is found in the previously cited co- 
pending U.S. Patent Appln. No. , entitled, "METHOD FOR RESIDUAL 
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FORM IN MOLECULAR MODELING," filed of even date. This computation fills in two 



partitions of the vector 



given previously. The first partition is called p , the kinematic 



\Pu) 

residual, and the second partition is called p u , the dynamic residual. The kinematic residual 
is computed from the difference between a q , which is passed-in from the (implicit) 
integration submodules 66, and the derivative computed by each joint: 

q-W(q)u = p q 

The dynamics residual is also computed. Starting with a given state of the 
molecular model, i.e., given (q,u,u) and the system parameters, a program routine models 
the 'environment' of the MBS. Such routines are readily available to, or can be created by, 
practitioners in the computer modeling field. The routine takes the values (q,u) determined 
by and passed in from the integration submodules 66 and returns (the state-dependent) 



T(k) = 



, the applied spatial force for a body k at its hinge point Q k , and cr(k) , the 



F(k) 

hinge torque for the body k. The dynamics residual, p u (k) , associated with generalized 
speeds u(k) for the body k is then computed by the following steps: 

1 . Perform the calculations for the molecular model by the Residual Form as described 
above with the passed-in state values (q,u,u) ; 

2. Generate f(k) , the spatial load balance for each body k in the model having n bodies: 



( N 



T{k) = M(k)A{k) + 



*> k {k)(le{k) N cD k {k)) 



-T(k) 



* = 1,. 

3. Compute p u (k) 



for k =n to 2 by - 1 

p u (k) = H(k)f(k)-a(k) 
I = inb(k) 

end 

p u {\) = H{\)f{\) 

The Residual Form method evaluates the extent to which the system 
differential equations are satisfied. Zero residual indicates that the applied forces are in 
balance with the inertia forces. However, this does not mean the system is in static 
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equilibrium, but rather that the applied forces would reproduce the given u when applied to 
the system in the state (g,u) . The residuals can be interpreted as that additional hinge torque 
needed to balance the applied and inertia forces. In the literature this method is known as 
either inverse dynamics, or the method of computed torques. It governs the case where the u 
are all prescribed. At this point all the computations required for the Residual Form are 
complete. The residuals p q and p u are used directly by the implicit integrator in the 
integrator submodule 68. 

Second Kinematics Calculations for the Molecular Model 

To carry out the Direct Form method, calculations in addition to the first 
kinematics calculations are required. These additional calculations are termed the second 
kinematics calculations. The values P(k),D(k), l y/ k (k), l K k (k) are computed as follows: 

1 . Perform the calculations for the Molecular Model by the Residual Form as 
described above, i.e., the first kinematics calculations. 

2. P(k) , the articulated body inertia of each body is initialized. 

P(k) = M(k), k = l 9 ... 9 n 

3. The objects below are then generated: 

for k = n to 2 by -1 

D(k) = H{k)P(k)H\k) 
G = P{k)H\k)D-\k) 
T = E 6 - GH(k) 
y k (k) ='/(*) r 
i K k (k)= i /(k)G 
i = inb(k) 

p(0 +=y k (k)P(k)Y\k) 

end 

D(l) = H(\)P(\)H\l) 
The functional dependence of these quantities is only upon the generalized coordinate q. 

Therefore, the first kinematics calculations are programmed in anticipation of performing the 

second kinematics calculations. 

Forward Dynamics Calculations 

Finally, u is calculated by sweeping inboard, then outboard, of the molecule: 
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z(k) = 0 6 , k = \,...n 
fork =n to 2 by-1 

s(k) = p u (k)-H(k)z(k) 

o(k) = D'\k)s(k) 

i = inb(k) 

z(i) +=V(fc)z(/:) + ^(^W 
end 

s(\) = p u (l)-H(l)z(l) 
u(l) = D- l (l)£(l) 
u(l) = u(l) 

for k = 2 to n 
i = inb{k) 

S{k) = '^'(kWO + HWuik) 
u(k) = u(k)- i K k \k)<5(i) 
end 

With the First and Second Kinematics Calculations, and the Forward Dynamics Calculations, 
the Direct Form method is available. 

Direct Form Method for the Equations of Motions 

The Direct Form method takes the current state (q,u) and computes the derivatives 
(q,u) using the above algorithms, which are then used by the integration method to advance 
time. 

Given: (q,u) 
Compute: {q,u) 

1 . Compute q using joint specific routine as above 

2. Perform above First Kinematics Calculations with u = 0 

3. Generate residuals p u as above 

4. Negate the residuals p u =-p u 

5. Perform Second Kinematics Calculations 

6. Compute it using Forward Dynamics Calculations above 

The Direct Form method produces the hinge accelerations u in response to the 
applied forces acting on the system. Fig. 5 summarizes the computation steps of the Residual 
Form method and the Direct Form method. 
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JACOBIANS IN IMPLICIT INTEGRATION 

The MD equations which model a molecule (such as a protein), are 
implemented as a multibody system (MBS). These equations represent Newton's Laws and 
are expressed as a set of differential equations y = f(y,t) . The differential equations are 
5 implemented using a suite of Order( N ) multibody dynamics methods. 

To advance the equations in time, in accordance with the present invention, an implicit 
method of numerical integration is used, in particular, L-stable implicit integration methods, 
such as implicit Euler, Radau5, and SDIRK3. 

An important ingredient of this integration process is formation of the 
10 Jacobian of the differential equations. This is 

Since the function /is itself computed by an algorithm rather than by an explicit formula, the 
Jacobian computation represents a substantial amount of work. In the simplest approach, the 
Jacobian can be formed numerically by differencing the derivative routine. This is a delicate 

1 5 operation because the quality of the Jacobian is a tradeoff between round-off and truncation 
errors. Typically half the working precision in the result is retained by choosing a good 
perturbation size in the difference scheme. In practice, though, this is difficult to do. 

However, the structure of the governing equations may be exploited to 
improve the Jacobian computation. The exemplary multibody dynamics methods illustrate 

20 this. The algorithms involved compute exact derivatives, even though numerical methods are 
used to execute the formulas. The derivatives obtained are in error by amounts that depend 
upon round off and the conditioning of the multibody system under consideration. But no 
approximations are involved at the equation level. 

In general, G, the iteration matrix used in the Newton loop of the implicit 

25 integrator has the form G = E - a J , where E is the identity matrix and a is be some scalar 
function of the timestep. See the previously referenced U.S. Patent Appln. No. 

, entitled "METHOD FOR LARGE TIMESTEPS IN MOLECULAR 

MODELING," filed of even date, for a description of implicit integrators. Changes in step 
size require refactoring G , but not reforming J . Reforming J is needed only when the 

30 Jacobian is needed at a new state. G is used in a linear subproblem within a Newton loop. 
The following is solved: 

Gky = -r{y% 
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where r(y) is the residual function for that particular implicit integration method. 

As shown later, J has a special structure, which is inherited by G . This 
means that equation solving with G can be done at reduced cost, if the structure is exploited. 

The quality of the Jacobian affects the ability to solve the nonlinear equations 
5 resulting from discretization in the integrator. Failure to solve the Newton loop may require 
retraction of a trail step and reduction of the integration time step. The timestep should be 
controlled by accuracy, rather than failures in the Newton loop. 

Computation of Analytic Jacobian 

The Jacobian J is a matrix which represents a linearization of the equations of 
10 motion. Normally, the governing equations for a dynamical system are linearized around an 
equilibrium state, or perhaps a state of steady motion. In this case, the equations are 
linearized around an arbitrary state so all possible contributing terms should be developed. It 
is customary to describe J in terms of its partitions: 

(dq dq} 

dq du 

dii dii 

K dq du j 

1 5 Structure of J qq AND J qu 

The q equation is q = Wu , where the matrix FT has a block-diagonal structure. 
Each block depends upon the joint type. Pins and slider joints give rise to 1 by 1 identity 
blocks. A ball joint generates a 4 by 3 block that expresses the Euler parameter derivatives in 
terms of Euler parameters and body angular velocity measure numbers (generalized speeds). 

20 From the q equation above, these two partitions of the Jacobian matrix are 

formed: 

dq 

J qu =w 

These equations are to be interpreted for symbolic purposes only. In practice, 
25 there is no need to generate the matrix ^explicitly. The non-zero block diagonal elements 
are filled in as discussed in the previous section on the kinematic residual. 



J(q,u) = 



J J ^ 

\J uq J uu J 
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Structure of J uq AND J uu 

d(-M' x p u ) 

The it derivatives are more complicated. Since u - -M p , — - 

8q 

d(-M~ x p\ 

and must be computed. (Note that p u is the residual of the dynamics routine 

du 

developed earlier). Here a key result from the field of Numerical Analysis is used to avoid 
the derivative of the matrix inverse. 

Suppose A(x)y(x) = b(x) , then we can write y(x) = A' l (x)b(x) . If 

dv( x} 

y(x 0 ) = z is known, and the value of — 1 x=x0 must be obtained, we have 

dx 



dy 
dx 



= A-\xO)-^-(b(x)-A(x)z) 

x=xO VX 



x=xO 

where z = A~ l b is held fixed when forming the right-hand side of the equation. 
10 Applied to the described multibody routines, 

^ = " £(^>„(<7>",0)) = -M~ x JLp u (q 9 u 9 z) 

where z = M~ x p u {q,u,0) . This result avoids the computing of the derivative of M~ l , which 

is a hypermatrix. The matrix inverse is "pulled-out". In the above equation, "x " is either the 
generalized coordinate q or the generalized speed w, depending on the Jacobian partition that 
15 is to be computed. 

In summary, to compute the blocks J uq and J uu , three steps are followed: 

1. Given (q,u) compute z = -M ~ l p u (q, w,0) using the Direct Form method. This 
simply says to compute u from the current state and save it as the variable z. 

2. Compute the analytic Jacobian of the dynamics residual routine. In this step, the 

20 matrix — p u (q,u,z) is to be formed. This quantity clearly depends upon the vector z 

ox 

computed in Step 1. Note that the numerical value of the residual p u in this step is 
zero for each element since we are computing the Jacobian around a consistent 
solution of the motion equations. The partitions J uq and J uu of the Residual Jacobian 

are obtained by substituting "q" and "u" for the "x" above. 
25 3. Back-solve the result of Step 2 with the mass-matrix to obtain the desired block. The 

back-solve operation is accomplished in the Direct Method routine by processing a 
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residual vector into a u vector. The Second Kinematics Step only needs to be 
performed once, since the back-solves are done at the nominal value of the state. In 
fact, the Second Kinematics routine must have been called in Step 1 while computing 
z, so the variables should still be cached. 
In words, the Jacobian of our derivative routine can be formed by back-solving the Jacobian 
of our residual routine. The residual Jacobian is much simpler to compute than the derivative 
Jacobian. Steps 2 and 3 above are derived in the following sections. 

The Residual Jacobian 

The computation of the Residual Jacobian is closely related to the Residual 
Form method for dynamics, which is summarized here: 

1 . Perform an outboard pass that computes the kinematic data that depends upon 
position and velocity. 

2. Make a call to the force routine which generates the atomic forces and consolidates 
them into spatial loads acting on the bodies. 

3. Perform another kinematic pass that computes acceleration level quantities (using 
passed-in u ), and combines inertia forces with the spatial loads from Step 2. 

4. Perform an inward pass that generates the residual at each joint. This pass recursively 
computes the resultant of the (spatial) inertial and applied forces for the 'nest' of 
bodies outboard of the joint in question. The residual is the projection of the resultant 
on the joint's degrees of freedom, given by the joint map data. 

At a high level the residual computation can be considered to depend upon 
two kinds of forces: 'motion forces' and external forces. The motion forces are computed 
directly by the multibody system. The external forces are available to the multibody system 
from a force modeling routine that computes the various interatomic forces such as 
electrostatics and solvents. A similar procedure is followed when computing the Jacobian. 
The multibody system builds the Jacobian of the motion forces, and combine it with the 
Jacobian of the external forces. 

Partitioning the Force Field into Effects 

There are many forces that may be acting on the molecule, and these forces 
may be computed in various intrinsic coordinate frames that are most convenient for that 
particular force "effect". For example, electrostatic terms may be computed using multipole 
methods and spherical coordinates, covalent terms may be computed in terms of torsion and 
bond angles, and solvent forces may be computed in global Cartesian coordinates. During 
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the computation of the Residuals, these forces are transformed from their intrinsic coordinate 
frame to the MBS coordinates. 

The same exchange occurs to compute Jacobians. The native Jacobians in 
their intrinsic coordinates are be brought into the MBS coordinates. This requires the use of 
5 the chain-rule to transform between intrinsic and the MBS generalized coordinates. It is 

important that each effect co-computes its function value and Jacobian, because many of the 
same terms are needed for each computation. Each effect is transformed into a set of spatial 
loads T eJTect (k) , where k is the index of a generic body in the system. The totality of these 

effects is given the symbol T(k) . 

10 Effect Jacobians Brought into the Residual Jacobian 

At a high level, the residual routine was previously implemented from the 

equation 

P u =H0(MA-T) 

The implementation uses Order( TV ) methods which are immediately obvious from the 
15 equation above. In this equation T is a vector of spatial loads acting on the pivots of the 
multibody system, where each element is a spatial load (a 6-vector composed of one force 
and one torque). It actually represents all effects other than inertia loads or pure hinge loads. 
The term in parentheses represents the load balance for each body. The first term is the 
inertia force, the next is the spatial load. M(k)A(k) is the spatial inertia force for a typical 
20 body. This is built from the body mass properties and the spatial acceleration of the body 
pivot. The spatial acceleration is computed before the residual routine is executed by the 
Forward dynamics routine. The operator H0 is implemented in a routine that performs an 
Order( TV ) inward pass. 

Even without knowing anything about the details of the computation implied 

dT 

25 by the equation above, the contribution of , the spatial load Jacobian (the effect Jacobian) 

dx 

to — - , the residual Jacobian, can be immediately inferred: 

dx 

dx dx 

The ... refers to terms not involving the effect Jacobian. Again, q or u for "jc", is 
substituted, depending upon which partition of the Jacobian is being computed. 
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The role of an effect Jacobian in the residual Jacobian is the same as the role 
of the effect in the multibody equations. This means that T contributes to the residual p u in 

the same way that a column of contributes to . Both are processed by the same 

dx dx 

operator H<f> . This is a crucial point because it means that no new method is needed for this 

dT 

part of the Jacobian computation. (A different method is need to obtain ). 



Thus, given the effect Jacobian, its contribution is assembled into the residual 
Jacobian by operating on it with the original residual routine, treating it as a multi-column set 
of spatial load vectors. This is a direct consequence of the linearity of the equations. 

The columns of the residual Jacobian play the same role in the derivative 
Jacobian routine as the residual vectors play in the Forward Dynamics routine. The dynamics 
routine performs a back-solve on a data vector it receives, and doesn't need to know what the 
data is, just what operation to perform on it. This applies to all the routines. 

Adding the . . . terms, the chain rule is used to show the whole equation: 



At this level, there are four contributions to the Jacobian: the inertia forces, the spatial forces, 
and contributions due to changes in the operators <fi and H. The quantities z and y refer to 

{MA - T) , and (f>z , respectively, which are held fixed while evaluating the last two terms. 

The numerical values of these terms are already available from the residual computation. 
Another observation about the above equation is that the operator <f> depends only upon q, but 

not u. Thus, this term remains constant while computing the partition J uu . Similarly, the 

spatial load can be split into its constituent effects, some of which do not depend upon w. In 
general, this means that knowledge of the multibody equations can be exploited to optimize 
the computations. 



described, but a description of how to form the term has not been made. These details are in 
the following sections. 

Computing the Effect Jacobians and Combining with the Residual Jacobian 

So far a high-level description of the Jacobian computation has been given. It 
can be seen that the computation has a very algorithmic flavor to it. There are very distinct 



dx 




Up to now, what to do with once the term has been computed has been 

dx 
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phases to the task, just as there were also for the Forward Dynamics routine described 
previously. There, computation of atomic forces is clearly the bottleneck step. Yet the 
overhead in the multibody equations for dealing with forces is fairly small. In that case, a 
call was made to the force routine, and what occurs inside the routine was ignored. When it 
5 comes to the Jacobian, this aspect is less true. 

A call is still made to obtain the Effect Jacobian, but there is a lot of 
processing needed before the Effect Jacobian can be assembled into the Residual Jacobian. 
The details of dealing with force fields to produce Jacobians are covered in the next sections 
and an example of incorporating electrostatics is developed. All other loads follow a similar 
10 development. 

Electrostatics as an Example Effect 

The basic premise of electrostatics is that the force between two charged 

particles is 



15 This is a classical inverse square law. F g , the force acting on particle i due to particle J 9 

depends upon the charge of the particles, attractive for oppositely charged particles, repulsive 
for like charges. The symbol r & is a unit vector directed from particle i to particle j; r i} is the 

distance between the particles; k is a unit-dependent constant related to the strength of the 

forces. Of course, the force acting on particle j is equal and opposite to that acting on the 
20 particle i. Hence, given a collection of particles interacting through Coulomb's Law given 

above, the net force acting on each particle is computed by summing the pair-wise forces. For 

this example, the forces are computed in global Cartesian coordinates. 

With the atomic forces in hand, the multibody forces can be generated. The 

system of forces acting on the particles of each body is replaced by a spatial load acting at the 
25 pivot of each body. The atomic forces are first expressed in a body- fixed basis, and then 

shifted to the pivot using the station coordinates of the particular atom to which the force is 

bound. 

F.. is a vector. The derivative of F & with respect to dr i} , small changes in the 
particle's relative positions is Q v , a tensor. It is such that dF i} = D^dr- . For Coulomb force 
30 the tensor is 
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Some observations: 

1 . The Force Jacobian is a matrix of size n atoms x n atoms . Each element is a 3 by 3 tensor. 

The (ij) block gives the derivative of force on atom i with respect to small changes in 
5 the position of atom j. In general, every force model is required to support an intrinsic 

Jacobian method for analytical processing. 

2. Storage requirements for the Force Jacobian quickly become impractical. This leads 
to the notion of interface "contraction" where the Jacobian of all the forces acting 
pair- wise between the atoms are reduced or "contracted" to acting pair-wise between 

10 the bodies. 

3. Because the Coulomb force is a pair-wise interaction, each force contributes to two 
blocks in the overall Jacobian. Thus, each force is processed at constant cost, and the 
overall Jacobian is computed at a cost proportional to the number of atoms squared, 
i.e., Order( N 2 ). This is the same as the computational cost of the force itself! This is 

15 a rather good result for computing analytic Jacobians. A numerical Jacobian requires 

a fresh force computation each time an element of the state is perturbed. This leads to 
cubic growth, i.e., Order( N 3 ), in the cost of the numerical Jacobian. Hence, the 
analytic Jacobian is much cheaper to compute as well as more accurate than a 
numerical Jacobian. 

20 4. The computation of the Jacobian is conveniently done in the same routine that 

computes the force (co-computation). However, it typically needs to be done far less 
often than the force computation. Therefore, a flag can be used to trigger the Jacobian 
calculation only when needed. 

Coupling to the Displacement Gradient 
25 Having obtained the (intrinsic) Force Jacobian, it is necessary to process the 

Jacobian further. This is due to the fact that the multibody system is formulated using 
relative coordinates. The chain rule is applied to each atomic force F(kJ) and is called 

"coupling to the displacement gradient". This denotes the global Cartesian force o\i atom i, 
which resides on body k. 

30 dFjkJ) = ^y dF(kJ) dr(p,s) 

tyfj £l~dr(p 9 s) dqj 

where r(p,s) is the position of atom "s" on body "p". 
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The first term in the sum selects an element of the Force Jacobian which was - 

just computed. The quantity ^ r (P> s ) - $ m e i ement Q f th e displacement gradient. A typical 

dq. 

term gives the change in an atom's position due to a small change in a generalized 
coordinate. Note that this term is strictly a kinematical quantity having nothing to do with the 
5 force computation. Thus, the Force Jacobian can be computed once and then continually 
reprocessed by the chain rule for each coordinate in the multibody system. This step 

dr. 

represents a matrix vector multiply, since — - is a column vector with n atona entries (each a 3 

vector), and the Jacobian is a square matrix n atoms x n atoms , where each element is a 3 by 3 
tensor. 

10 It is possible to improve this computation, since many of the entries in the 

displacement Jacobian are known to be zero. This is due to the fact that incrementing a 
particular hinge does not displace every atom in the system, but only those outboard of the 
displaced hinge, and not disjoint from the hinge in question. For instance, rotating the base 
body induces a change in all atoms of the system. But perturbing a torsion angle on any 

15 terminal body induces a change only to those atoms resident in the terminal body. Therefore, 
roughly speaking, about half the work may be saved by optimizing the computation. This 
reduction comes from a strictly knowledge-based approach. 

Interface Contraction 

In the process of forming T(k) , the spatial load on body k 9 the load comes from the atomic 
20 forces acting on atoms that belong to body k. Each force is transformed from global to local 
coordinates, and then shifted to the body pivot. A concise statement of this procedure is: 

The operator (f>{k^i) is 



~£ 3 


p(k,i) 


' N C k (k) 


2? 


_% 


13 




N C k {k) 



0(k,i) = 

25 where p(k,i) are the fixed station coordinates of atom i on body k. Note that the new 

quantity T(k,i) appears. This is just the atomic force turned into a spatial load at the atomic 
position of atom i: 



T(k,i) = 



F(k,i) 
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The first element of the atomic spatial load is zero because there is no torque exerted by the 
force field on individual atoms. 

Now, T(k) relates atomic forces to body spatial loads. So, the derivative of 
this equation relates differential atomic forces to differential spatial loads: 

j til d( lj d<lj 

= Tt(k) + T 2 (k) 

The second term, T 2 (k) , in this equation is discussed at the end of this section, as it involves 
the spatial loads, but not the load derivatives. This means the term can be treated generically, 
without worrying how the spatial loads were computed. 
1 0 Substituting the definition of T(k) , into T^k) : 

ie* p=l sep 

dr(p,s) dqj 



nbod 



dT{k) dr(p,s) 



tl?pdr(p,s) dqj 

where now the symbol zZiQ £ V^(& >t ) d?XM) - g m e i ement Q f t jj e row-reduced Force 

or(p,s) JS dr(p,s) 

Jacobian. This Jacobian relates differential (body) spatial loads directly to differential atomic 
displacements. Again, r(p,s) refers to the global position of the atom s on the body p. The 

1 5 term— is formed by summing each atomic Force Jacobian element into the destination 

dr(p,s) 

element in the reduced Jacobian, weighted by the atoms' 0(kJ) matrix. Each element of the 
row-reduced Jacobian is a 6 by 3 matrix. Hence the rows of the Force Jacobian have been 
contracted. The contraction is evident in the notation: the numerator has only a body index, 
while the denominator has both a body and an atom index. Depending on the number of 
20 atoms per body, the row reduction can provide a savings in both storage and execution time 
when differential spatial loads must be formed. 

Note that the row-reduction procedure only needs to be done once before 
computing the Residual Jacobian. The overhead of performing the reduction is more than 
offset by the reduced cost of the smaller matrix vector products which must be formed. 

25 Note that in forming d7X&) t k ere j s no neec j to save ^ e i emen t s Q f t h e atomic Force 

dr(p 9 s) 

Jacobian. That is, each element ^H^ill on iy ne eds to be available while its contribution to 

dr(p y s) 
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r ; \ is being computed. So, more than one element of the big Force Jacobian is not 
dr(p,s) 

required at a time. 

The global coordinates of a typical atom, r(p 9 s) , are computed in terms of 
r(p) , the global coordinates of body p's pivot, and p(p,s) , the station coordinates of the 
5 atom: 

r(p,s) = r(p) + N C k (p)p(p 9 s) 
By differentiating we find, (with some results from the kinematics of finite rotations): 

Augmenting this equation with the additional equation X(p,s) = X(p) and defining w, the 
10 spatial derivative 

" X 

w = 



20 



8r 
Idq. 



the result is: 





X(p,s) 




' Ap) ' 


w(p,s) = 


dr(p,s) 




dr{p) 




L *j J 




L ^ J 



15 The vector X can be interpreted as generating a rate of change of orientation for each body. 
It is a field quantity, in the sense that it can potentially vary at each point in space. For rigid 
bodies undergoing pure rotations (without deformation), it is constant for each body affected 
by the rotation. An Order( N ) algorithm for computing A recursively is described in a latter 
section. 

With the application of the above equation to the equation for T x (Jc) , the final 
reduction formula is reached: 

nbod 



p(dw(p) dq. 

8T(k) ±y\®3 6T(k) 
dw(p) ^|_0, dr(p,s) 
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Each element of — ^-^ , the reduced Jacobian, is now a 6 by 6 matrix. An element of the 
dw(p) 

reduced Jacobian relates the spatial load at the pivot of a body to a spatial derivative 
occurring at another pivot in the system (after coupling to the displacement gradient of the 
pivots). 

5 The row reduction consolidated all the atomic forces on each body, leaving the 

spatial load derivatives coupled to displacement derivatives of all the atoms in the system. 
The column reduction consolidated all the atomic displacement derivatives, leaving the 
spatial load derivatives coupled to spatial pivot derivatives. Thus, the Force Jacobian is 
recast into a "native" form. Working with the reduced Jacobian speeds up computation of the 

10 spatial load derivatives by roughly the square of the number of atoms per body. This speedup 
can easily approach a factor of 100 or more. 

This section shows, using electrostatics as an example, how to build an 
atomic-level Force Jacobian. This Jacobian relates differential atomic forces to differential 
atomic displacements. With the introduction of the idea of interface contraction, differential 

15 spatial loads are shown to be related to differential atomic displacements through a row- 
reduced Force Jacobian. Another improvement to the computation finally results in a fully 
contracted Jacobian that relates differential spatial loads to differential spatial displacements. 

Now, the Force Jacobian must be coupled to the spatial displacement 

dr(p) 

gradient described above. A matrix vector multiply performs this step, and scales with the 
20 number of bodies in the system, not the number of atoms. 
The second term, T 2 (k) , is given by: 

iek dqj 



and involves the original spatial loads T(k) and derivative of the operator 0(k,i) 





P(k,i)\ 







o, 



C*(k) 



25 Thus 



T 2 (k)= X 



N dC k (k) p(k,i) N dC k (k) 
N dC k (k) 



O3 



T(k,i) 



where 



N dC k (k) = N dC k (J) 'C k (k)+"C k (i) l dC k (k) 
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is computed recursively from the base body outward and 



'dC k (k) = ^f-^-dg(k) k = \,...n 



where dq(k) is defined in the next section, and 



d'C'jk) 



, the partial derivative of the 



interbody direction cosine matrix is a function of the joint type connecting the bodies: 

d'C k (k) 



pin 
slider : 



d <lk 
d'C k (k) 



= -K sin(q k ) + A cos(q k ) + AJL sin(q k ) 



= 9* 



ball and free : 



d'C*(k) 
ds x 



= 2 



0 








-2s x 


~ £ 4 


*3 




-2e x 



ds, 



d l C k (k) 



= 2 



= 2 



—2s, 



-s A 



'-2s* 



-2s, 



d l C k (k) _ 



ds 4 



= 2 



-2s* 



0 -s x 
s t 0 



In the next Section the Force Jacobians are combined with the Inertia Force Jacobians to 
finally form the Jacobian of the Residual Routine. 

The Residual Jacobian 

The previous Section describes the formation of the Force Jacobian — - , 

dw(p) 

10 which must be coupled to the spatial displacement gradient, in order to form the derivative of 
the spatial forces. This section describes the formation of the spatial displacement gradient 
and the formation of the Jacobian of the residual routine. 

The recursive algorithms for computing the entire Jacobian are described. The 
Jacobian algorithm is not actually set up to compute the Jacobian. As is typical of automatic 

15 differentiation routines, it computes the matrix vector product J uq dq + J uu du for arbitrary 

passed-in values of the vectors dq and du. In practice, to compute the Jacobian, the "Jacobian 
Routine" is effectively called repeatedly with a series of Boolean vectors (a vector with one 
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entry set to 1 and all other entries set to zero.) Each call generates the corresponding column 
of the Jacobian. Note that some of the steps have already been or are being computed for the 
Residual Form method or the Direct Form method (the Forward Dynamics Calculations), but 
are reproduced here for clarity. 

1. Given (q,u) compute z = -M ~ l p u (q,u 9 0) using the Direct Form method. Also set 

ii — z and recompute A{k) , then p u , which recomputes f(k) . 

2. Perform contraction to compute the fully row- and column-reduced Force Jacobian, 
dT{k) 



dw(p) 



as described in section, "Interface Contraction' 



dT ^ dT . 
— = o — <t> 
dw dr 



10 Steps 3 through 10 below are used to fill the columns of : 

3. Compute the analytic Jacobian partitions of the q terms: 

qq dq 
J = W 

using joint routines similar to those needed for the First Kinematics Calculations. 

15 4. Compute q derivatives of position quantities and terms for spatial gradient: 

Previously described methods were used to fill in certain joint-specific fields. These 
quantities consisted of l C k (k) , the interbody direction cosine matrices, r QiQk (k) , the 
spanning vector for each body, and H{k) , the joint map for each body's inboard joint. 
To refer to the derivative of each of these quantities, a prefix 'd' is added to the 

20 symbol name to make this reference generically. Thus, l dC k {K) means the derivative 

of the direction cosine matrix l C k (k) . 

Each interbody direction cosine matrix (and all the joint- specific) quantities depend 
only on the generalized coordinates of an individual joint. Thus, l dC k (k) is nonzero 
only when the derivative is taken with respect to any of the coordinates for body k. 
25 To properly 'seed' the recursions being studying, a vector dq is passed in to the 

routine. For Jacobian computation we set one entry is set to 1 , and all the other 
entries to 0. Then, the needed preliminary quantities are generated by a typical loop: 

i dC k (k)= d '^ k) dq(k) k = l,...n 
dq(k) 
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The partial derivatives of the direction cosine matrices are generated analytically and 
displayed in the section, "Interface Contraction" above. These partial derivatives do 
not depend upon the particular column of the Jacobian that is being computed. 
Setting a particular entry of dq to 1, and all the rest to 0 has the effect of annihilating 
the correct subset of the seed quantities. 

, the partial derivative of the interbody spanning vector is given by 

dq k 

dr QtQk {k) 



dr QlQ > (k) 

dq k 
dr Q ' e '(k) 

dq k 



dr QiQ >(k) 



= A>(k)> slider 
= g x4 , ball 

= 0j, pin 

0 0 0 0 1 0 0 
0 1 0 



free 



0 0 0 0 
0 0 0 0 0 0 1 
A(k) here refers to the body's sliding axis which connects it to its parent body. 

dH(k) 



-, the partial derivative of the joint map is 



dH(k) 
dH(k) 

dH(k) 
d <lk 



= Qs> pin, slider 
= [0, 0,], ball 

d'c k (k) 



free 



With the above definitions of the partial derivatives, the recursions are seeded with 
the following loop: 

for k = 1 to nbod 
i dC k (k) = d ' Ct(k) dq(k) 

dr Q - Q >(k) = drQ ' 6>(k) dq(k) 

dH(k)=^dq(k) 
end 
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After execution of these loops, all bodies have l dC k {k) , dr QiQk (k) , and dH(k) , their 
interbody derivative quantities available. 

One new quantity needed in the spatial displacement gradient computation is 
also computed. This is A(k) from the section on Interface Contraction, the rotation 
axis that generates the rate of change of orientation for each body outboard of a 
perturbed joint. Here, this variable is given the symbol d0(k) , the differential 
rotation axis for each body is 

for k = 1 to nbod 

d0{k) = A(k) dq{k\ pin 

d0(k) = O 3 , slider 

d0(k) = not needed for ball, free 

end 

Since arbitrary perturbations to a set of Euler parameters do not produce a pure 
rotation, column contraction cannot be used when computing the corresponding 
column of the Jacobian. The row-reduced Force Jacobian can still (and must) be 
used.. 

After seeding the recursions, N dC k (k) 9 dr OQk (k), l d<fr k {k) y dA(k) are computed: 

N dC k (l) = i dC k (l) 
dr OQ >(\) = dr Q > Q '{\) 

dA(l) = d0(l) 
for k - 2 to nbod 

N dC k (k)= N dC k (iYC k (k)+ N C k (iYdC k (k) 
dr OQk (k) = dr OQk (z) + N dC k (i)r Q G> (k) + N C k (i)dr QiQ * (k) 
'a b 
c d 

dA(k) = 'C**(*)cU(0 + d0(k) 

end 
where 

a = l dC k (*), b = df Q & (k) 'C* (k) + r aa l dC k (k), 
c = g 9 d = 'dC'ik) 

Compute q derivatives of velocity: 

A loop that computes the rate of change of joint velocity due to change in joint angle 
starts the process: 



l d^ k {k)^ 
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for k = 1 to nbod 
i dv k (k) = dH\k)u(k) 
end 

This quantity is the rate of change of joint velocity due to change in joint angle. 

Obviously, it is nonzero only for joints whose map contains coordinate dependence. 

For free joints, the generalized speeds produce relative linear velocity that depends 

5 upon the joint orientation. 

After computing i dv k (k) , dV{k) , the derivative of the spatial velocity of each 

body, is computed. This is done by the following loop: 

dV(\) = i dv k (\) 
for k — 2 to nbod 

dV(k) = l d<f> k *V(i) + WdV(i) + i dv k {k) 
end 

6. Couple Force Jacobian to spatial displacement gradient to compute T x (k) 

for k = 1 to nbod 

f^dw(p) dqj 

end 

7. Compute second term of the Force Jacobian T 2 (k) and append to T x (k) : 

fork= 1 to nbod 

" N rfC* (A:) ^(A:, I) "</C* 
O3 N dC k (k) 



dT(k) = W)*^ 

IE* 

end 



8. Compute q derivatives of acceleration-related terms: 

Again the process starts with a loop that computes l da k (k) = dH\k)ii{k) : 
for k = 1 to rtfeorf 
15 i da k (k) = dH\k)u(k) 

end 

This is the change in joint acceleration due to a change in coordinate. Then, dA(k) , 
the derivative of the spatial acceleration of each body is computed. 

dA(l) = i da k (l) 
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'V *(*)-» 



V(k) 



V(A:) 

N a>\kj 
% a (*) 



, 'dV k (k) 



dV(k) 



'da> k (k) 
'dv e >(k) 

N da> k (k) 
N dv Q >(k) 



where — > means the 6 vector is decomposed into two 3 vectors, 
for k = 2 to nbod 

dt\ = 'dC k * (k) ( N co k (i) x ( N a> k (i) x r aa (£))) + 
\ N da> k (0 x ( N co k (0 x r aa (A:))) + 
'C**(A:) ( N a>" (0 x ( N dot (0 x r aa (£))) + 
J W (i) x ( W (0 x </r aa 

cfo = 'd0 k '(k)A(i) + *0 k '(k)dA(i) + jj^j 

^'^'C**(A:)^*(0 
</*>/ = WC** (*) W (i) + 'C** (*) (0 
A 2 = dcojx V (&) + ay x '</<»* (A:) 
rfr3 = 2rf*y j x f v a (A:) + 2a>j x Wv a (A:) 
"^2" 



dA(k) = da+ ' + 'da k (k) 
\_dt3 

end 

The symbols introduced here with = are meant to be temporary variables not needed 
after computation of dA(k) . 

After computing the spatial acceleration derivatives, the computation of 

df(k) , the spatial inertia force derivatives, is performed: 

for k = 1 to nbod 

dv\ = N da> k (k) x ^ (A) ■ "at + N a> k (k) x ^ (A) ■ N da k 

dv2 = N dco k (k) x ( N co k (k) x p(k)) + N Q) k (k) x ( N deo k (Ar) x p(A:)) 



dT(k) = M(k)dA(k) + 
end 



dvl 
dv2 



dT(k) 



Compute dp(k) , the joint residual derivative for body k: 
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for k = nbod to 1 

dp(k) = dH(k)f(k) + H(k)df(k)-der(k) 
i = inb(k) 
if i >0 

df(i) = df(i) + 'd0 k (k)f(k) + '/ (k)df(k) 

end 
end 

After executing this routine, the values stored in dp(k) are the new column of the 

3 

Residual Jacobian — p u (q, u,z) . 

dq 

5 10. Back-solve the — result of previous step with the mass-matrix to obtain the desired 

dq 

du. 

dq dq 

The back-solve operation is accomplished in the Direct Form method routine by 
processing a residual vector into a u vector. The Second Kinematics Calculations 
10 only needs to be performed once for the whole Jacobian, since the back-solves are 

done at the nominal value of the state. In fact, the Second Kinematics routine must 
have been called in Step 1 while computing z, so the variables should still be cached. 
Steps 1 1 through 13 below are used to fill the columns of J uu : 

1 1 . Compute u derivatives of velocity: 

15 This routine takes a passed-in vector du and computes i dv k {k) = H\k)du(k) . Then, 

dV(k) , the derivative of the spatial velocity of each body, is computed: 

dV(l)= i dv k (l) 

for k = 2 to nbod 

dV(k) = l <p k \k)dV{i) + l dv k (k) 

end 

12. Compute the velocity-induced derivative df(k) . As presented here, the routine is 
specialized for the case of no velocity dependent external loads. The surviving terms 

20 are those due to changes in inertia forces alone. Even if there were changes in 
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external loads, it would only be required to include the contribution of dT(k) as 
before. 



dA{\) = 



O3 
0, 



for k = 2 to nbod 



dtl = 'C k '(k) 



'( N do> k (i) x( N a> k (1) x r Q & (*))) 
( N a> k (i)x( N dco k (i)xr Q < Q >{k)))j 

da = '0 k '(k)dA(i) + 



O3 
dtl 



d&>j= i C k '(k) N da> k (i) 
dtl = dcojx i <o k (k) + a>jx 'da> k (k) 
dt3 = Ida? j x ' v a (k) + lcojx 'dv e ' (k) 
dtl 

dA(k) = da + 



dt3 



end 



After computing the spatial acceleration derivatives, df(k) , the spatial inertia force 
derivatives, is computed: 

for k = 1 to nbod 

dvl = N da> k (k) x (k) ■ N a? k + N a? k (k) x (k) ■ N dco k 
dvl ± N d<o k (k) x ( "at (k) x p(k)) + N a>" (k) x( N da> k (k) x p(k)) 

dvl 



dT(k) = M(k)dA(k) + 



end 



dvl 



13. Compute dp(k) , the joint residual derivative for body k: 

for k = nbod to 1 
dp(k) = H(k)df(k) - da(k) 
i = inb(Jk) 
ifi>0 

df(i) = df(i) + */ (k)df{k) 

end 
end 
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After executing this routine the values stored in dp(k) are the new column of the 

d 

Residual Jacobian — /?„(^«,z). 
du 

14. Back- solve the — result of previous step with the mass-matrix to obtain the desired 
du 

du. 
du' 

5 ®L = _ M - X dp 

du du 

The back-solve operation is accomplished in the Direct Method routine by processing 
a residual vector into a u vector. The Second Kinematics Calculations only need to 
be performed once, since the back-solves are done at the nominal value of the state. 
In fact, the Second Kinematics routine must have been called in Step 1 while 
10 computing z, so the variables should still be cached. 

The above steps complete the computation of the analytic Jacobian as long as 
the forces only have dependence on q . This accommodates the classical situation where all 
atomic forces are derivable from a potential function. To accommodate velocity-dependent 
forces, such as simple viscous damping, some of the above steps need to be modified as 
15 follows: 

In Step 2 above, we also need to compute the contracted velocity Jacobian 

— — - , which is block diagonal, must also be computed. 
dr(Jk) 

In Step 6 above, the computation of T x (k) must be augmented with the 

contracted velocity Jacobian: 

for k = 1 to nbod 

end 



where 



dr(k,i) N ,„ k 



= N dC k (k)\ N v Qk (k) + N <v k (k) x p(k,i)~] + 
dqj J 

N C(k)[ N dv Q > (k) + N da> k (k) x p(k, I)] 
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A step is then added after Step 11, which is called Step 11a. This new step 
computes dT(k) : 

l^dr(kj) duj 



where 



duj 



= N C k (k) [ N dv Qk (k) + *W (*) x p(k,i)] 



While executing Step 12 above, the last loop for df(k) is modified by 
subtracting the velocity-dependent force derivative dT(k) : 

for k = 1 to nbod 

dvl = *W (k) x ^ (A:) • N co k + "V (A:) x (k) • 
rfv2 = N dco k (k) x("a> k (k)x p(k)) + N a> k (k) x( N da> k (k) x p(Ar)) 

"rfvl" 



^r(Ar) = M(k)dA(k) + 



</v2 



-dT(k) 



end 

The rest of the Steps remain the same. 
10 Fig. 6 summarizes the operational steps of the Analytic Jacobian method, 

which has been described in detail above. 

Fig. 7 shows a plot of the accuracy of the numerical Jacobian versus the 
accuracy of the analytic Jacobian for an exemplary MD system. In the best case in which the 
perturbation was perfectly selected, the digits of accuracy for the generalized coordinates (q) 
15 and generalized speeds (w) from the numerical Jacobian, illustrated by line 152, were still 
only half that of the analytic Jacobian, illustrated by line 150. 

ADDITIONAL EMBODIMENTS 

The present invention has many embodiments besides the examples described 
above. The list below has other embodiments and applications: 

20 • Order of Forces included in Jacobian 

Any order of the forces to be included in the Jacobian, include, but not limited 
to Order( N ), Order( N 2 ), Order( N 3 ), and Order( N 4 ). An example of an Order( N ) force 
field would be an electrostatic force field using fast multi-pole expansion methods (see, for 
example, Greengaard, The Rapid Evaluation of Potential Fields in Particle Systems, 

25 Massachusetts Institute of Technology Dissertation, 1988) rather than direct method which is 
Order(W 2 ). 
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• Analytic Jacobian for Direct Form 

When the governing equations are in Direct Form, the so-called "forward 
dynamics" form of the equations is obtained. In this form, the equations process a state 
vector and applied efforts and generate the acceleration at each of the joints modeled in the 
5 system. 

u = A/"' (/) 

The Jacobian then represents the partial derivatives of the accelerations with respect to 
elements of the state vector. The preferred embodiment shows several algorithmic methods 
for computation of these partial derivatives. The methods are exact and do not utilize 
10 numerical approximations to form derivatives. 

The Direct Form produces the u partitions of the Jacobian: 

*(*-'(/)) 



Juq ' dq 



and 



a(M- (/)) 

du 

15 by using an algorithmic counterpart to the function which computes the u function. 

The computation of M' x f is accomplished in Order( N ) floating point 
operations (flops) by utilizing the operational inverse of the Mass Matrix: 

M' x ={I~ H^K] D' 1 [I - H^Kf 

where all the block matrices are defined previously in the Second Kinematics Calculations, 
20 and each factor represents an operator that can be applied to an w-vector in Order( N ) flops. 

Since u = [/ - H'VK) ZT 1 [/ - H^Kf f , from the chain rule: 

J uq = [/ - H*¥K]D- 1 [I - HVKf df/dq + 
d ([/ - HWK]) I dq D~ l [I - HVKf f + 
[I-HVK] d(D- l )/dq[l-fP¥K] T f + 
[I-HVK]D- 1 d([l -H"VKY)ldq f 

and similarly for J uu . 

Thus the present invention can be used for equations of the Direct Form in 
25 producing algorithms that compute each of the above terms in closed form. 
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Hence the present invention provides many advantages. Analytical Jacobians 
are much more accurate (with twice the number of significant digits). Computing from the 
Residual Form instead of Direct Form is much more efficient. The "Contraction" of rows 
and columns from "number of atoms" to "number of bodies" reduces the size of the force 
5 Jacobian matrices. Jacobian computations are of the same order as computation of forces, 
rather than an extra order higher if each column has to be perturbed. Thus, if the forces are 
computed in Order( N 3 ) operations, for example, a numerical Jacobian requires Order(7V 4 ), 
whereas an analytic Jacobian requires only Order( N 3 ) operations. By controlling the range 
of loop structure in Jacobian calculations, computations can reduced even further (just 
10 compute for outboard bodies). 

Therefore, while the foregoing is a complete description of the embodiments 
of the invention, it should be evident that various modifications, alternatives and equivalents 
may be made and used. Accordingly, the above description should not be taken as limiting 
the scope of the invention which is defined by the metes and bounds of the appended claims. 
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WHAT IS CLAIMED IS: 
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1 . A method of modeling the behavior of a molecule, comprising 
selecting a torsion angle, rigid multibody model for said molecule, said model 

having equations of motion; 

selecting an implicit integrator; and 

generating an analytic Jacobian for said implicit integrator to integrate said 
equations of motion so as to obtain calculations of said behavior of said molecule. . 

2. The method of claim 1 wherein said analytic Jacobian is derived from 
an analytic Jacobian of the Residual Form of the equations of motion. 

3. The method of claim 2 wherein said analytic Jacobian J comprises 



(dg 


dg 


dg 


du 


du 


du 


K d 1 


du 



11 J<p* 
. J J 



;and 



j egjJWu) d± = w 

dg dg qu du 

j u =?»=-M-> d ^ q > U > z) and J uu = *L = - M -> dp ^ q ' U > z) 
dq dq du du 

where q are the generalized coordinates, u are the generalized speeds, Wis a joint map matrix 
1 5 and M is the mass matrix and p u is the dynamic residual of the equations of motion, and z is 
-M-'p u (q,u,0). 



4. The method of claim 3 wherein said implicit integrator selecting step 
comprises an L-stable integrator. 

5. A method of simulating the behavior of a physical system, comprising 
20 modeling said physical system with a torsion angle, rigid multibody model, 

said model having equations of motion; and 
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integrating said equations of motion with an implicit integrator; said implicit 
integrator having an analytic Jacobian to obtain calculations of said behavior of said physical 
system. 

6. The method of claim 5 wherein said analytic Jacobian is derived from 
5 an analytic Jacobian of the Residual Form of the equations of motion. 

7. The method of claim 6 wherein said analytic Jacobian J comprises 



dq 


dq] 




dq 


du 


( 


du 


du 


V 


dq 


du j 




dq 


d{Wu) 



Jqq Jqu 



; and 



and J = %- = W 



qv dq dq qu du 

J m =*=-M-* &&&*> and j m -*-- M - 1 

q dq dq du du 

10 where q are the generalized coordinates, u are the generalized speed, W is a 

joint map matrix and Mis the mass matrix and p u is the dynamic residual of the equations of 
motion, and z is -M~ x p u (q 9 u, 0) . 

8. The method of claim 7 wherein said implicit integrator comprises an 
L-stable integrator. 

15 9. Computer code for simulating the behavior of a molecule, said code 

comprising 

a first module for a torsion angle, rigid multibody model of said molecule, said 
model having equations of motion; and 

a second module for an implicit integrator to integrate said equations of 
20 motion with an analytic Jacobian to obtain calculations of said behavior of said molecule. 

10. The computer code of claim 9 wherein said analytic Jacobian is 
derived from an analytic Jacobian of the Residual Form of the equations of motion. 
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1 1 . The computer code of claim 10 wherein said analytic Jacobian J 

comprises 



J = 



r dq_ 
dq 
du 



{dq 




; and 



J. 



dq 
dq 



d(Wu) 
dq 



and J, 



gu 



dq 



W 



J, 



uq 



du 
"dq 



-i dp u (q 9 u,z) 
dq 




du 
du 



-i dp u (q 9 u 9 z) 
du 



where q are the generalized coordinates, u are the generalized speed, FT is a 



joint map matrix and M is the mass matrix and/> w is the dynamic residual of the equations of 
motion, and z is -M~ l p u (q,u,0) . 

12. The computer code of claim 1 1 wherein said implicit integrator 
comprises an L-stable integrator. 

13. Computer code for simulating the behavior of a physical system, said 
code comprising 

a first module for a torsion angle, rigid multibody model of said system, said 
model having equations of motion; and 

a second module for an implicit integrator to integrate said equations of 
motion with an analytic Jacobian to obtain calculations of said behavior of said system. 

14. The computer code of claim 13 wherein said analytic Jacobian is 
derived from an analytic Jacobian of the Residual Form of the equations of motion. 

15. The computer code of claim 14 wherein said analytic Jacobian J 

comprises 
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J = 



f dq_ d£\ 

dq du 

du du 

ydq du j 



r J J \ 



; and 



dq = d(Wu) 
qq ~ dq~ dq 



and J = — = fV 
du 



j a «,-jr' te^ and 

dq dq du du 

where # are the generalized coordinates, m are the generalized speed, is a 
joint map matrix and Mis the mass matrix andp w is the dynamic residual of the equations of 
motion, and z is -M~ l p u (q 9 u 9 0) . 

16. The computer code of claim 15 wherein said implicit integrator 
comprises an L-stable integrator. 

17. A method of modeling the behavior of a molecule, comprising 
selecting a torsion angle, rigid multibody model for said molecule; and 
generating a partition of an analytic Jacobian J for said model, said partition 

selected from the group comprising J gq9 J qU9 and J UU9 and 



'dq dq} 
dq du 
dti du 



K dq du j 

where q are generalized coordinates, q are derivatives of said generalized 
coordinates with respect to time, u are generalized speeds, and u are derivatives of said 
generalized speeds with respect to time. 



18. The method of claim 17 wherein said model has equations of motion 



and 
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J m = * m -M-> 2aSS&*l and J m = *=-M-**** 9 ^' 
m dq dq du du 

where Wis a joint map matrix defining relationships between q and u 9 Mis 
the system mass matrix defining relationships between u and the generalized forces on said 
bodies, p u is the dynamic residual of said equations of motion, and z is -Af ~ l p u (q 9 u 9 0) . 

19. The method of claim 1 8 wherein said molecule comprises a number of 
atoms and wherein said selecting step further comprises 

reducing said number of atoms to a smaller number of rigid bodies so as to 
reduce the number of computations for Jacobian partitions J uq and J uu . 

20. The method of claim 19 wherein said number of computations is 
reduced approximately by the average number of atoms in each of said rigid bodies. 

21. A method of modeling the behavior of a physical system, comprising 
selecting a torsion angle, rigid multibody model for said physical; and 
generating a partition of an analytic Jacobian J tor said physical model, said 

partition selected from the group comprising J qq , J qU9 J uq and J uu , and 



dq 


dq 


dq 


du 


du 


du 


dq 


~du~ 



r J J \ 

\^ uq ^ uu J 



where q are generalized coordinates, q are derivatives of said generalized 
coordinates with respect to time, u are generalized speeds, and u are derivatives of said 
generalized speeds with respect to time. 



motion and 



22. The method of claim 21 wherein said physical model has equations of 



J 



dq d(Wu) 
dq dq 



and J =^- = W 

qu du 



dq dq 



du 



du 
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where Wis a joint map matrix defining relationships between q and u, Mis 
the system mass matrix defining relationships between u and the generalized forces on said 
bodies, p u is the dynamic residual of said equations of motion, and z is -M~ x p u (q, w,0) . 



23. Computer code for simulating the behavior of a molecule, said code 

comprising 

a first module for a torsion angle, rigid multibody model for said molecule; 
and 

a second module for generating a partition of an analytic Jacobian J for said 
model, said partition selected from the group comprising J qq9 J qU9 and J UU9 and 



(dq 




dq 


du 


da 


du 


K dq 


du j 



uq J uu J 



where q are generalized coordinates, q are derivatives of said generalized 
coordinates with respect to time, u are generalized speeds, and u are derivatives of said 
generalized speeds with respect to time. 



motion, and 



24. The computer code of claim 23 wherein said model has equations of 



dq dq qu du 



j =^L = -M^ d ^' U ' z) and J uu =^L = -M^ dp ^ z) 
dq dq du du 

where Wis a joint map matrix defining relationships between q and w, Mis 
the system mass matrix defining relationships between it and the generalized forces on said 
bodies, p u is the dynamic residual of said equations of motion, and z is -M~ l p u (q 9 u 9 0) . 

25. The computer code of claim 24 wherein said molecule comprises a 
number of atoms and wherein said first module further comprises 

code reducing said number of atoms to a smaller number of rigid bodies so as 
to reduce the number of computations for Jacobian partitions J uq and J uu . 
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26. The computer code of claim 25 wherein said number of computations 
is reduced approximately by the average number of atoms in each of said rigid bodies. 



27. Computer code for simulating the behavior of a physical system, said 
code comprising 

a first module for a torsion angle, rigid multibody model for said physical 

system; and 

a second module for generating a partition of an analytic Jacobian J for said 
model, said partition selected from the group comprising J qq , J gU9 J uq and J uu , and 



J = 



dq 


dq 


dq 


du 


du 


dti 


dq 


du 



1 ' Jqq Jqu 



where q are generalized coordinates, q are derivatives of said generalized 
coordinates with respect to time, u are generalized speeds, and u are derivatives of said 
generalized speeds with respect to time. 



motion, and 



28. The computer code of claim 27 wherein said model has equations of 



j A JJw«) a 

w dq dq qu du 



j - dii - M -i d Pu(<l>">z) md j _ { dp u (q,u 9 z) 

dq dq du du 

where Wis a joint map matrix defining relationships between q and w, A/is 
the system mass matrix defining relationships between u and the generalized forces on said 
bodies, p u is the dynamic residual of said equations of motion, and z is -M~ l p u (q,u 9 Q) 



29. The method of claim 1 further comprising calculating said analytical 
Jacobian in no more than order N p * ] computations, where N is proportional to the number of 
rigid bodies in said model, and wherein a calculation of a numerical Jacobian for said model 
requires order N p computations. 
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30. The method of claim 5 further comprising calculating said analytical 
Jacobian in no more than order N p l computations, where N is proportional to the number of 
rigid bodies in said model, and wherein a calculation of a numerical Jacobian for said model 
requires order N p computations. 

3 1 . The computer code of claim 9 wherein said second module further 
comprises code for calculating said analytical Jacobian in no more than order N**" 1 
computations, where N is proportional to the number of rigid bodies in said model; and 
wherein a calculation of a numerical Jacobian for said model requires order N p computations. 

32. The computer code of claim 13 wherein said second module further 
comprises code for calculating said analytical Jacobian in no more than order 
computations, where N is proportional to the number of rigid bodies in said model; and 
wherein a calculation of a numerical Jacobian for said model requires order N p computations. 

33. The method of claim 1 wherein a calculation of a numerical Jacobian 
for said model has an accuracy of K digits and further comprising 

calculating said analytical Jacobian for said model with an accuracy of at least 

2K digits. 

34. The method of claim 5 wherein a calculation of a numerical Jacobian 
for said model has an accuracy of K digits and further comprising 

calculating said analytical Jacobian for said model with an accuracy of at least 

2K digits. 

35. The computer code of claim 9 wherein a calculation of a numerical 
Jacobian for said model has an accuracy of K digits and said second module further 
comprises 

code for calculating said analytical Jacobian for said model with an accuracy 
of at least 2K digits. 

36. The computer code of claim 13 wherein a calculation of a numerical 
Jacobian for said model has an accuracy of K digits and said second module further 
comprises 

code for calculating said analytical Jacobian for said model with an accuracy 
of at least 2K digits. 
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37. The method of claim 1 wherein said analytic Jacob ian is derived from 
an analytic Jacobian of the Direct Form of the equations of motion. 

38. The method of claim 37 wherein said analytic Jacobian J comprises 



J = 



'dq d£\ 

dq du 

du du 

K dq du t 



; and 



dq d(Wu) 
Bq dq 



and J o =^- = W 

qu du 



= a(M-(/)) j(M->( f) ) 

dq uu du 



where q are the generalized coordinates, u are the generalized speeds, W is a 
joint map matrix defining relationships between q and u 9 M* 1 is the operational inverse of the 
system mass matrix M which defines relationships between u and the generalized forces / on 
said bodies. 

39. The method of claim 5 wherein said analytic Jacobian is derived from 
an analytic Jacobian of the Direct Form of the equations of motion. 

40. The method of claim 39 wherein said analytic Jacobian J comprises 



J = 



' dq_ dq^ 

dq du 

du du 

^dq du j 



qq - qu 



\ J uq 



; and 



uu J 



dq d(Wu) 
qq ~dq~ dq 



and J =^- = W 
qu du 



= a(M-(/)) 

dq du 
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where q are the generalized coordinates, u are the generalized speeds, Wis a 
joint map matrix defining relationships between q and u, M" 1 is the operational inverse of the 
system mass matrix M which defines relationships between u and the generalized forces / on 
said bodies. 

41 . The computer code of claim 9 wherein said analytic Jacobian is 
derived from an analytic Jacobian of the Direct Form of the equations of motion. 



42. The computer code of claim 41 wherein said analytic Jacobian J 



comprises 



J = 
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dq) 
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du 
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du 
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du j 
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( Jqq Jqu 
\^ uq ^ uu J 



; and 



dq dq 



and J qu =^- = W 

gu du 



j(M->(f)) j(M-(f)) 
dq uu du 



where q are the generalized coordinates, u are the generalized speeds, W is a 
joint map matrix defining relationships between q and u, M* 1 is the operational inverse of the 
system mass matrix M which defines relationships between u and the generalized forces / on 
said bodies. 

43. The computer code of claim 13 wherein said analytic Jacobian is 
derived from an analytic Jacobian of the Direct Form of the equations of motion. 



44. The computer code of claim 43 wherein said analytic Jacobian J 



compnses 
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qu 



; and 
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dq dq du 

J -~ dq andJ -- £ 

where # are the generalized coordinates, u are the generalized speeds, Wis a 
joint map matrix defining relationships between q and w, M" 1 is the operational inverse of the 
system mass matrix M which defines relationships between u and the generalized forces / on 
said bodies. 
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METHOD FOR ANALYTIC JACOBIAN COMPUTATION IN MOLECULAR 

MODELING 

ABSTRACT OF THE DISCLOSURE 
A method for obtaining analytic Jacobians used in implicit integration 
methods in the computations for the dynamics of a physical system. With this method, the 
Jacobian with at least twice the number of digits of accuracy as a numerical Jacobian can be 
computed. This also results in the implicit integration method being more efficient because a 
smaller number of iterations are required to solve the nonlinear stage equations of the 
equations of motion, as well as the ability to take larger timesteps. This speedup in 
computation is very useful in molecular modeling. 
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